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<3\ ' Abstract. The only way to map the Galaxy's gravitational potential 

^(x) and the distribution of matter that produces it is by modelling 
£h ' the dynamics of stars and gas. Observations of the kinematics of gas 

I ^ . provide key information about gradients of within the plane, but lit- 

tle information about the structure of out of the plane. Traditional 
Galaxy models assume, for each of the Galaxy's components, arbitrary 
flattenings, which together with the components' relative masses yield 

t— I ■ the model's equipotentials. However, the Galaxy's isopotential surfaces 

should be determined directly from the motions of stars that move far 

^? ■ from the plane. Moreover, from the kinematics of samples of such stars 

\ that have well defined selection criteria, one should be able not only to 

t-h ■ map <& at all positions, but to determine the distribution function /j(x, v) 

C of each stellar population i studied. These distribution functions will con- 

tain a wealth of information relevant to the formation and evolution of 
the Galaxy. An approach to fitting a wide class of dynamical models to 
■ the very heterogeneous body of available data is described and illustrated. 

Oh; 

6 ■ 
a: 

1. Introduction 

Models of the Milky Way can be usefully divided into three classes: (i) pho- 
tometric models, (ii) kinematic models, and (hi) dynamical models. Models of 
the first type were pioneered by Kapteyn and Shapley, a more recent example 
is given by Bahcall & Soneira (1980). These models specify the distribution of 
stars without saying anything about their motions. The models of Bienayme 
et al. (1987) and Ratnatunga et al. (1989) are typical kinematic models: they 
specify both the distribution and the motions of stars. However, in these mod- 
els the velocity distribution at each point must be somehow specified without 
fully exploiting Newton's laws of motion. In practice the Galaxy is decomposed 
into 'components' such as the thin disk, the thick disk, the metal-poor halo etc, 
and at each point in space each component is assigned an ellipsoidal, usually 
Gaussian, velocity distribution. One might hope to use the Jeans' equations to 
couple the ellipsoids assigned at neighboring points, but Eddington (1915) al- 
ready showed that this is feasible only if the potential is of Stackel's special form, 



1 Invited talk presented at the meeting "Formation of the Galactic Halo ... Inside and Out", 
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which it almost certainly is not. So the velocity ellipsoid at each point in each 
component must be arbitrarily assigned at risk being dynamically impossible. 
In consequence, models of this type may have a role to play as representations 
of observational data, but they by no means fully exploit the potential of stellar 
surveys. 

This contribution emphasizes the value of constructing models in which 
stellar dynamics plays an integral role, and explains why such models have until 
recently been of an extremely limited nature, and why the time is now ripe for 
them to come into their own. 

2. Orbits and Integrals 

The basic result of galactic dynamics is that galaxies can be considered to be 
made up of orbits in a smooth gravitational potential. From this it follows that 
a prerequisite for dynamical modelling is an understanding of such orbits. In 
the simplest circumstances orbits in a steady three-dimensional potential are 
characterized by three 'isolating integrals,' that is functions ij(x, v) of position 
and velocity that are constant along any orbit. Through the work of Lindblad 
(1933), Contopoulos (1963), and Henon & Heiles (1964) it has long been clear 
that most orbits in typical axisymmetric galactic potentials admit three such 
integrals. Two integrals are available from elementary mechanics, I\ = E, the 
star's energy, and I2 = L Z , the star's angular momentum about the potential's 
symmetry axis. The problem is that, for most systems of interest, no general 
formula exists for a third integral, 1%. By Jeans' theorem, the simplest assump- 
tion one can make on an equilibrium system is that the probability density 
/(x, v), that a given star will be found at the phase-space point (x, v), is a 
function f(Ii, l2,Iz) of three isolating integrals. Hence ignorance of I3 makes 
it impossible to express a general distribution function (DF) / in its simplest 
form. 

External galaxies have been modelled dynamically more extensively than 
has the Milky Way by the simple stratagem of ignoring 1% (e.g., Binney et 
al. 1990). When the observational data are sparse, this procedure can lead to 
acceptable models. Unfortunately, it has been known since Jeans' classic (1915) 
paper that this stratagem cannot lead to an acceptable model of the Milky 
Way, because it predicts that the velocity dispersions in the radial and vertical 
directions must be equal; near the Sun they differ by nearly a factor of 2. 

2.1. The Oort Lindblad approximation 

Oort (1932) and Lindblad (1933) decomposed the stellar motion into a vertical 
oscillation and the motion parallel to the plane. The latter can be further 
divided into azimuthal rotation and radial libration. Let (R, if, z) be cylindrical 
polar coordinates, then the difference between a star's energy E = \v 2 + $(x) 
and the energy E c + Er associated with motion parallel to the plane provides 
an approximate form of ^3 : 

E Z = E-(E C + E R ) = \v\ + *{R, z) - *(R, 0). (1) 
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Here E C (L Z ) is the energy of a circular orbit with angular momentum L z and 

L 2 

E R = lv 2 R + ^ eS (R)-E c (L z ) where $ eS (R) = $(J2, 0) + ^- (2) 

is the energy associated with radial librations. This Oort-Lindblad form of I3 
plays a key role in all published determinations of the local mass density within 
the disk (e.g., Bahcall 1984, Bienayme et al. 1987, Kuijken & Gilmore 1991), 
and readily explains the Schwarzschild velocity ellipsoid (e.g., §4.2 of Binney & 
Tremaine 1987). In fact it enables one to explain the significant skewness of the 
observed distributions of v v velocities (Cuddeford & Binney 1994). 

However, (|l|) is only an approximation and its limitations give rise to the 
leading uncertainty in the column density of the disk near the Sun. The problem 
is that when E z is used as third integral, one necessarily finds that (vrv z ) = 
with the result that even away from the plane the velocity ellipsoid is aligned 
with the R- and z-axes. By contrast, if the Galaxy's potential were spherically 
symmetric, the velocity ellipsoid would everywhere align with the direction to the 
galactic center. We can be sure that the real situation is intermediate between 
these two extremes, leading to a fundamental uncertainty of ~10%-20% in the 
local surface density (Kuijken & Gilmore 1989). 

This inability of the O-L approximation to yield an unambiguous determi- 
nation of the local vertical force K z (z) and thus the disk's surface density is the 
more disappointing as the latter is one of the very few parameters of the galactic 
potential that is primarily determined from the O-L approximation. Indeed the 
O-L approximation is applicable only to stars on nearly circular orbits - disk 
stars. These probe the potential only near the plane, and for most purposes 
the interstellar medium (ISM) provides a more readily observed probe of the 
potential near the plane, since its radio-frequency spectral lines can be observed 
throughout the disk without hindrance from the dust which heavily obscures 
even moderately distant disk stars. Consequently, the Galaxy's circular-speed 
curve v c (R) has long been determined from observations of the ISM, and these 
data can in principle be used to determine K z near the plane (Merrifield 1993, 
Malhotra 1994). 

The only way to improve substantially on the knowledge of the galactic 
potential that we have gleaned from observations of the ISM is to study the 
dynamics of stars that move far from the plane and thus are high- velocity stars 
when in the solar neighborhood. The O-L approximation does not apply to 
these stars, and some other approach to ^3 must be developed. 

2.2. Approaches to I 3 

How can we obtain adequate models of the orbits of high- velocity stars? One pos- 
sibility is directly to integrate the equations of motion. For each initial condition 
chosen, this produces a stream of phase-space positions [x(t^), v(ifc)], k = 1, N. 
Two tasks must be addressed before a useful galaxy model can be based on this 
raw material: (i) characterize the orbits in some systematic way, such as assign- 
ing values of E,L Z ,I%; (ii) put the output data onto some sort of grid so that 
one can subsequently decide whether a given orbit will eventually pass through 
a given point x, and, if so, with what velocity v. 
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A recent paper by Zhao (1996) exemplifies this approach. Zhao adopts a 
simple model of the rotating potential of the Milky Way's bulge and calculates 
several hundred orbits in it. He characterizes his orbits by energy and the time- 
averages along them of the angular momentum about the potential's long and 
spin axes. He uses the orbit's positions and velocities to determine associated 
occupation probabilities and velocities in each of 1000 spatial cells. 

Merritt & Fridman (1995) employ a very similar technique to model ellip- 
tical galaxies. The main difference between their work and Zhao's is that they 
follow Schwarzschild (1993) in characterizing orbits not by time-averaged angu- 
lar momenta, but through their initial conditions. In this technique all orbits of 
a given energy are launched from one or more two-dimensional surfaces. These 
surfaces are carefully chosen so that (a) any orbit can be obtained by launching 
from some point of one of them, and (b) as far as possible, each point on the 
surfaces generates a distinct orbit. I2 and I3 are simply the values taken by two 
convenient coordinates for the starting surface(s) at the orbit's launching point. 

In Oxford we have developed an entirely different technique which involves 
thinking of orbits as three-dimensional surfaces in six-dimensional phase space 
rather than as time-series (Kaasalainen & Binney 1994, and references therein). 
These surfaces are topologically equivalent to tori, so we call this the "torus 
method". At present the only three-dimensional orbits we can construct are 
axisymmetric ones. In this case each orbit is characterized by three special 
isolating integrals, the actions J r ,J[,J v . Here J r is the generalization of the 
energy Er of equation (|2|); J; is the generalization of E z in equation (||); and 
J^p = L Z is simply the angular momentum about the symmetry axis. 

In place of the time-series (x/%, v/%) above we obtain x, v as Fourier series in 
the so-called angle variables 6 r ,6i,Q(p. Along an orbit angle variables have the 
remarkable property of increasing linearly in time: 6i{t) = 9i(0) + 0,{t, where the 
J7j are constants characteristic of the torus. Consequently, the probability that 
a star will be observed at any point on its torus is uniform in the variables 

The actions are special in the following important sense. The phase-space 
volume occupied by orbits with actions in an elementary cube in integral space 
of size <i 3 J is / d 3 x<i 3 v = (27r) 3 d 3 J. In other words, when one uses actions as 
integrals, equal volumes in integral (or 'action') space correspond to equal vol- 
umes in the full six-dimensional phase space. Since the distribution function 
/(J) is a probability density in phase space, (27r)~ 3 / (J) is by this result also 
the probability density in action space: fi(J)d 3 J is the probability that a star 
of population i has actions in <i 3 J. This fact makes it easy to understand the 
physical implications of any particular form of the DF fi of each population. 
Consequently, when actions are used as integrals one knows from general astro- 
physical considerations in advance of detailed modelling roughly how fi should 
depend on its arguments. In fact, we will specify the functional form of fi and 
fit the data by merely adjusting a small number of parameters in the specified 
form. 

3. Dynamical Models of High-Velocity Populations 

Several high-velocity stellar populations may be distinguished spectroscopically, 
for example, RR-Lyrae stars, blue- horizontal-branch stars, and low-metallicity 
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subdwarfs. A dynamical model would associate to each of these sub-components 
a total number of stars N{ and a DF fi (I\ ,h,h) normalized such that / d 3 x d 3 vfi 
=1. Equipped with a model of this type, one could replicate within the com- 
puter any well-defined observational selection, i.e. predict the survey probability 
distribution for any survey. 

For example, Flynn & Freeman (1993) observed the radial velocities of a 
sample of stars that were photometrically selected to be M giants located about 
20kpc from the Sun towards the south galactic pole. One can evaluate the 
model's corresponding survey probability distribution by selecting stars from 
the model according to precisely the same selection criteria and determining the 
distribution of radial velocities of these stars, and compare it with the observed 
distribution. The model would also predict the distribution of the stars' proper 
motions, and these could be compared with any measured values. Notice that 
in this comparison uncertainties in the distances to individual stars, that arose 
from, for example, a broad distribution in absolute magnitudes amongst the se- 
lected stars, would not play an important role so long as the absolute-magnitude 
distribution had been correctly modelled. 

It is easy to see that a model of the type just described unambiguously 
predicts the survey probability distributions for a bewilderingly large number of 
observational surveys. Indeed, in whatever direction and in whatever magnitude 
range one selects stars of a given photometric or spectroscopic type, both the 
radial- velocity and proper-motion distributions are determined once <!>, Ni and 
fi have been chosen. 

Of course, not all observations can be used to test the model; some will be 
employed to determine Ni and fi and there will be nothing to be learnt by 
using the model to replicate these particular observations. So it is important to 
estimate how many observations will be needed to establish a model. 

The model depends upon a handful of three-dimensional objects: 3>(x), and 
the fi(Ii, I-ii Iz)- If we assume (for the moment) that the Galaxy is axisymmetric, 
then $ is specified by the ISM throughout the plane and we have to guess only 
how $ changes away from the plane. For any such guess we can determine fi 
in a large part of three-dimensional integral space merely by determining the 
velocity distribution of each population i near the Sun (e.g., May & Binney 
1986). If the velocity distribution can be determined at one or two other well- 
chosen points in the Galaxy, we can complete our knowledge of fi. So it should 
be possible to predict the outcome of infinitely many very different surveys from 
the results of just a handful of surveys! Of course, our predictions will only be 
as good as our guessed form of <£. But it is clear that we have very much less 
freedom in extending this function away from the plane than we potentially have 
information from even a moderate number of surveys. 

The reason why even limited observational work leads to more information 
than the models have freedom is that Jeans' theorem renders the Galaxy a three- 
dimensional object, but one which we are privileged to observe in six-dimensional 
phase space. Even when our observations of phase space are seriously deficient, 
so that, for example, we determine neither the radial velocities nor the distances 
of stars, we still observe the Galaxy in more than three dimensions. Conse- 
quently dynamics renders our observations highly degenerate, and they become 
stringent tests of our model. 
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The simplest, yet rather complicated, dynamical models for the stellar halo 
will start by assuming that the Milky Way is both axisymmetric and fully mixed. 
The first assumption is certainly incorrect for the centre of the Galaxy, but this 
might have only little influence at large galactocentric distances. Also, as we 
have learnt at this workshop (Majewski, this volume), the assumption of a well 
mixed halo is an over-simplification. However, the size of both effects, non- 
axisymmetry and non-mixedness, is unknown and can most easily be measured 
by comparing the data to such idealised models. 



4. A Simple Example 

To demonstrate the viability of our Oxford approach and to illuminate the con- 
struction of our models and the prediction of survey probability distributions, 
we created a simple toy model and computed its predictions for a NGP proper 
motion survey. 

4.1. The Gravitational Potential 

In order to simplify the procedure we have chosen a scale free mass distribution 
consisting of a spheroid and a disk 

' nB -'- ] -1.4xl0 5 (i? 2 + [z/0.8] 2 )-^ 2 + 5.8xl0 5 i?^sech 2 — ^— (3) 



M pc- 3 v 0.06R 

with 7 = 1.8 and R, z measured in parsec. The spheroid's ellipticity is E2, 
while the disk's vertical structure is given by the usual sech 2 model with scale 
height z<i = 0.03-R. The parameters are chosen such that at Rq = 8kpc we have 
u c = 200kms -1 , A' 2 (l.lkpc)/(2vrG) = 74M pc" 2 , and z d = 240pc in agreement 
with observations. The local mass density is 0.013 and 0.055 M pc~ 3 for 
spheroid and disk, respectively. The gravitational potential is evaluated by 
inserting the ansatz (with polar coordinates r, #) 



$(r,0) = 4vrGr 2 " 7 



2088^% In cosh ^ 
pc 6 0.06 



(4) 



into Poisson's equation and solving for g($) using a multipole expansion. The 
ansatz (0) is chosen such that already a moderate number of multipoles gives 
the potential and forces to high accuracy. 

The advantage of choosing a scale invariant mass model is the existence of 
scaling relations, which reduce the number of independent dimensions of action 
space to two. 

4.2. A Distribution Function for the Metal-Poor Halo 

It is well known that for the most gravitational potentials surfaces of constant 
energy in action space are nearly planes A = aJ r + bJ[ + J v with a and b being 
functions of E. In the scale-free case, moreover, the shape of the energy surface 
is invariant under the scale transformation, hence a and b are constant. For our 
mass model we find a ~ 1.66, 6~ 1.16 and A ex £K 4_ t)/( 2 [ 2 ~">']) _ 
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An approximate DF depending on energy is obtained with f(A). If / ex. A* 1 , 
then the system has radial density profile poc^ 3 / 2-1- ^ 4-7 )/^ 2-7 ^. From this 
consideration we finally choose our stellar-halo DF to be 

/ halo (J)=^-( 6 - 7 )/( 4 - 7 )[vl + L circ (lkpc)] 2 ^/( 4 -^, (5) 

which yields pocr -7 for r<Clkpc and pcxr~P for r^>lkpc, where we have taken 
,3 = 3.5. 

Replacing A in / by A' = a'J r + b' ' Ji + J v produces a flattened and/or aniso- 
tropic model by shifting stars over energy surfaces. For our example we choose 
a' = 6' = 1.2, which mainly shifts the stars to slightly higher J r , i.e. to radially 
extended orbits. 

4.3. A Distribution Function for the Stellar Disk 

Unlike the mass model, our model for the distribution of population I stars will 
be an exponential disk. Suppose all disk stars were on exactly circular orbits. 
Then the disk's DF would be of the simple form 

/(J) = fo(Jtp) S(J r ) $(Ji), (6) 

and the mass dM with angular momentum between J v and X, + dJ v would be 
(27r) 3 /o( J<p) dJcp. For an exponential disk, this must equal MdjR\ exp(— R/Rd)R 
dR. Expressing R by the radius of the circular orbit with angular momentum 
Jip, one thus finds 

, Mj (dlt\ ( RcUJ r ) \ m 

We can get a vertically and radially warm disk by replacing the product of delta 
functions in (^) by (Binney 1987) 

0.(J)ty(J) ( J r^r(J) JA(J)\ , fi x 

2—2 — exp 2 2 — y°) 

afaf \ a$ of J 

with the orbital frequencies Slj(J). The functional form of (Ti(Ju>) determines the 
scale height as function of radius. Using the observed relation af oc exp(— R/Rd) 
and a similar one for a r , we finally get 



/disk(J) = (2n)m^(0)af(0) nr{m(J) 



x exp 



RtiTc(J(p J 



Rd 



j r n r (j) JA(j) 

a 2 (0) + af(0) 




(9) 



We have chosen the of (0) to result in <7 u = 36kms -1 and a w = 19kms _1 at Rq. 
Note that this simple form for the DF of a dynamically warm, exponential disk 
is valid for any axisymmetric potential. 
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Figure 1. Contours of the probability for observing proper motions 
Ho (in GC direction) and /U90 (in rotation direction) as computed 
from our model for a NGP survey with apparent magnitude limit 
14 < B < 18. The distributions for disk (left) and halo (right) are shown 
separately. The contours are spaced by 0.4 dex. 



4.4. An Example for a Survey Probability Distribution 

A survey probability distribution (hereafter SPD) resulting from our dynamical 
model is easily evaluated by a Monte Carlo technique as follows, (i) Choose 
(x, v) randomly from the ff, (ii) assign stellar age, mass, and metallicity accord- 
ing to some simple models for the SFR, IMF, and t-Z relation; (iii) use stellar 
evolutionary models to find absolute magnitude and colours; and (iv) 'observe' 
the star by computing the apparent observables and applying the selection cri- 
teria of the survey to compare with. Alternatively to using theoretical models 
of stellar evolution one could use some colour-magnitude relation as observed, 
say, for star clusters to assign 'stars' to the phase space points. 

As an example we computed from the above model for halo and disk the 
SPD of proper motions in the NGP direction. The SFR was taken to set in 12 and 
14 Gyr ago for disk and halo, respectively, with exponential decay times of 12 and 
2 Gyr. The common IMF was assumed to be of the form m~^ x+1 ^ where x = 1.35 
above and £ = 0.3 below O.8M with lower cutoff at 0.15M Q . We assumed 
a gaussian Z-distribution whose mean grows exponentially in time from 10 -4 
initially to 0.04 today. Stellar colours and magnitudes were otained using the 
models of Rezini & Voli (1981) and Maeder (1981 & 1991). To translate (x, v) 
into distance and proper motion we used zq = 7pc and v & = (9, 12, 7)kms _1 . 

The resulting SPD for the disk (Fig. |], left) clearly shows the effects of 
asymmetric drift: the steady downward shift in the centres of contours with 
increasing dispersion (contour size). This effect becomes much more significant 
for more local surveys (not shown). Though non-rotating, the halo SPD (Fig. 
|l[ right) peaks near zero proper motion, since most halo stars are too far away 
to show significant proper motion; the few nearby halo stars create the tail at 
H9o < 0. 
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Mo [mas/yr] fi [mas/yr] 



Figure 2. As Fig. [I] but for the joint distribution (solid) of disk and 
halo with a local normalization lOOphalo = Pdisk- Left and right panel 
are for 14 < B < 18 and 18 < B < 22, respectively. The dashed contours 
show the contribution due to the disk alone. 

In the joint distribution (Fig. left) the halo stars contribute negligibly. 
Thus, in order to use the halo stars as dynamical tracers, one has either to go 
for deeper surveys (Fig. [2], right), or to disentangle disk and halo stars - for 
instance in our model there are hardly disk stars in the colour-magnitude range 
B>12+10(B-V). 

5. Conclusions 

Very considerable advantages flow from modelling the Milky Way in a way that 
exploits Jeans' powerful theorem. It is especially useful to model stellar popu- 
lations that move far from the plane. What has held up this type of modelling 
for decades has been the difficulty of handling orbits that are constrained by a 
non-classical integral such as 1%. Several viable approaches to this problem are 
now available so that it should soon be possible to compare the new data that 
will flow from exciting observational tools such as Hipparcos and the 2dF with 
fully dynamical models. 

We have described some preliminary models based on the Oxford torus 
technique. In this method each population is assigned a DF /j that has a simple 
functional form, which describes the population's dynamical characteristics in 
astrophysically comprehensible terms. The parameters in fi are adjusted so as 
to optimize the fit between real catalogues and simulated ones. In principle any 
catalogue that has well-defined selection criteria can be simulated. 

The modelling process starts by guessing the Milky Way's potential &(R, z). 
Significant errors in this guess will lead to discrepancies between the real and 
simulated catalogues. Since the fitting process is strongly over-determined, we 
are confident that we will be able to determine $ to reasonable accuracy by 
iteratively adjusting it until the real and simulated catalogues agree. 
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